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Abstract. Delivering the full benefits of first principles calculations to battery 
materials demands the development of accurate and computationally-efficient 
electronic structure methods that incorporate the effects of the electrolyte 
environment and electrode potential. Realistic electrochemical interfaces 
containing polar surfaces are beyond the regime of validity of existing continuum 
solvation theories developed for molecules, due to the presence of significantly 
stronger electric fields. We present an ab initio theory of the nonlinear dielectric 
and ionic response of solvent environments within the framework of joint density- 
functional theory, with precisely the same optimizablc parameters as conventional 
polarizablc continuum models. We demonstrate that the resulting nonlinear 
theory agrees with the standard linear models for organic molecules and metallic 
surfaces under typical operating conditions. However, we find that the saturation 
effects in the rotational response of polar solvent molecules, inherent to our 
nonlinear theory, are crucial for a qualitatively correct description of the ionic 
surfaces typical of the solid electrolyte interface. 
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1. Introduction 

The development of better batteries is a critical step towards reducing energy reliance 
from oil to renewable energy resources. Experimental studies have historically 
made great progress in identifying promising battery materials pQ, but a number 
of challenges remain. For instance, the solid electrolyte interfaces (SEI) [5] of the 
electrodes play a crucial role in the thermodynamics and kinetics of battery operation, 
but these surfaces are currently poorly understood. Experiments have yet to even 
conclusively determine the compositions of many of these surfaces [3] . Other important 
challenges include understanding reaction mechanisms at the surfaces of electrodes and 
identifying new electrode materials. These problems are well-suited for computational 
study, which can complement an experimental approach through inexpensive, rapid 
but accurate calculations. 

However, electrochemical systems pose a unique challenge for theoretical studies: 
the processes of interest occur at an interface that requires simultaneous quantum- 
mechanical and statistical treatment. The electrode and reactants on its surface need 
to be described with a quantum-mechanical method in order to capture the level of 
detail required to predict chemical reactions. The liquid electrolyte plays an equally 
important role in determining the reaction pathways, and necessitates a statistical 
treatment due to the need to sample the configuration space of the liquid. 

The most straightforward approach to a combined statistical and quantum 
mechanical calculation is ab initio molecular dynamics [J], which is expensive since 
adequate statistical sampling necessitates several thousands of steps at the electronic 
structure level of detail. This cost may be ameliorated by combining classical 
molecular dynamics with electronic structure only for relevant parts of the system, as 
in the Quantum Mechanics / Molecular Mechanics (QM/MM) methods [5]. However, 
statistical sampling issues and the need for coupling constant integration for estimating 
free energies complicate the analysis of the results of any molecular dynamics based 
method. 

The complexity and computational cost due to statistical sampling can be avoided 
by working directly with equilibrium properties of the system. Joint density- functional 
theory (JDFT) [5] is an exact variational principle for the free energy of an electronic 
system in contact with a liquid, in terms of the densities of the two subsystems. This 
framework enables systematic approximations such as combining electronic density- 
functional theory for the system of interest with classical density-functional theory for 
the liquid environment. 

Polarizable continuum models (PCM's) [7] are a class of highly efficient simplified 
theories where the effect of the fluid is captured by placing the electronic system in an 
appropriately chosen dielectric cavity, optionally with corrections for physical effects 
such as cavitation energies and dispersion interactions. However, the efficiency of these 
models comes at the cost of empiricism and a loss of key physical features of the fluid. 

The empiricism of PCM approaches has been partially mitigated by constructing 
variants of the model [SJ |5] that are highly simplified approximations within the 
framework of JDFT. So far, PCM approximations 6, 8, 9, 10 have replaced the fluid 
with a linear dielectric response which turns out to be adequate for the solvation of 
most molecules and some surfaces, such as those of metals. However, the highly polar 
surfaces typical of battery systems impose strong electric fields on solvents that invoke 
a highly nonlinear response; linear response approximations lead to qualitatively 
incorrect results as we demonstrate in Section [3.41 
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In this paper, we present a systematic framework (Section |2.1[ ) for developing 
PCM-like approximations within joint density-functional theory, and use it to 
construct a nonlinear polarizable continuum model (Sections |2 .3| and 2.4) that is both 
inexpensive and sufficiently accurate to account for complex reactions, including those 
occurring on ionic surfaces. We show that the nonlinear dielectric model reproduces 
molecule solvation energies (Section 3.2), and with the inclusion of nonlinear ions, 
potentials of zero charge for metallic surfaces (Section 3.3) with accuracy similar to 
that of the linear model. Finally, we demonstrate that the inclusion of nonlinear 
dielectric saturation effects facilitates accurate predictions for ionic surfaces in solution 
(Section 3.4), making this model particularly suited for theoretical studies of battery 
materials. 



2. Nonlinear polarizable continuum model 

2.1. Joint density- functional theory framework for polarizable continuum models 

The fundamental quantity of interest for ab initio studies of electrochemical and other 
solvated systems is the free energy of a quantum mechanical system in equilibrium 
with a liquid environment. Therefore, the most direct route to this quantity is a 
theory in terms of the equilibrium densities of the two subsystems. Joint density- 
functional theory (JDFT) [Q\ is based on an exact variational principle for this free 
energy in terms of these equilibrium densities, and provides a rigorous framework for 
the development of practical approximations. 

The total free energy of such a system may be exactly partitioned as 

AjdftM-ZVJ] = Ahk N + $i q [{ N a }] + AA[n,{N a }] , (I) 

electronic liquid coupling 

where Ahk is the exact Hohenberg-Kohn electronic density functional [11], $i q is the 
exact free energy functional of the liquid [T3], and the remainder, A A, is the free 
energy for the interaction of the two systems. Minimizing the above functional yields 
the ground state electron density n(r) and the set of nuclear densities {N a (r)} for 
the fluid. 

In practice, each of the in-principle exact pieces of ([I]) needs to be approximated, 
and the power of the framework lies in the capability of independently selecting 
the level of approximation for each piece depending on the type of system, desired 
accuracy and available computational resources. The electronic system may be treated 
within the Kohn-Sham formalism [T3] with any of the standard exchange-correlation 
functionals, or if necessary, with correlated quantum chemistry methods or quantum 
Monte Carlo methods as demonstrated in [14"] . 

The liquid free energy may be treated within the rigid molecule classical density- 
functional theory formalism |15j . with an approximation for the excess free energy of 
the liquid; reliable functionals for liquid water have been constructed from its equation 
of state [13 [E] and functionals for other liquids are available in the literature. (See 
|17] for a survey.) The interaction of the two subsystems, A.A, may be treated using 
a density-only electronic density functional approach |18j . These approximations may 
be independently improved or simplified, as required for the system of interest. 

Using a classical density-functional theory for the liquid within JDFT is a powerful 
tool for studying solvated electronic systems. However, the complexity of the theory 
can occasionally obscure an intuitive physical interpretation of the results. This 
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intuition may be better obtained from simpler and possibly less accurate versions 
of the theory that capture the bare minimum of physical effects required to describe 
the systems and properties of interest. 

Polarizable continuum models (PCM) are highly simplified theories that account 
for liquid effects by embedding the electronic system in a dielectric cavity. The 
linear response approximation in PCM, however, is inadequate for the study of 
electrochemical systems that involve liquids in strong electric fields. Here, we develop 
a general framework for constructing PCM-like approximations within joint density- 
functional theory, which we use in the following sections to construct a nonlinear PCM 
with the same optimizable parameters as those of the linear model. 

We start by dividing the liquid contributions to the free energy functional into 
physical effects assumed to be separable in polarizable continuum models, and rewrite 
the last two terms of ([I]) in the following form[|] 

A dic i = $!q + AA 

= A e [s,e] + A K [s,fi] +JdrJ dr'^^ (p el (r) + + A cav [s}. (2) 

Dielectric response dominates the electrostatic interaction of a fluid consisting of 
neutral molecules alone, and the first term A e captures the corresponding internal 
energy. In addition to neutral solvent molecules, electrolytes typically include charged 
ions that contribute an additional monopole response. The optional term, A K , 
accounts for the internal energy of the ions if present in the solution. The densities 
of the molecules and ions of the solvent are modulated by the cavity shape function 
s(r), which in turn is determined by the electron density n(r). 

The third term of ([2| is the mean field electrostatic interaction of the liquid 
bound charge density pi q with itself and the electronic system of total charge density 
p c \. Here, p\ q = p e +p K includes dielectric and ionic contributions, while p c \ = n-\-p nuc 
includes contributions from the electrons and the nuclei (or pseudopotential cores) of 
the subsystem treated using electronic density-functional theory. The contributions 
from all remaining effects of the fluid, such as cavitation and dispersion, are gathered 
into the final term of ^ , A cav , and are assumed to depend only on the shape of the 
cavity s(r). We detail specific approximations for each of these terms in the following 
subsections. 

So far, the dielectric and ionic responses are still fully general, except for the mean- 
field assumption in their interaction with each other and the electronic system. In 
reality, these responses are nonlinear as well as nonlocal, while conventional polarizable 
continuum models [HI 13 [HI El EH assume both linearity and locality. In this work, we 
retain the local response approximation, but develop a nonlinear theory for dielectric 



and ionic response in sections 2.3 and 2.4 We obtain a linear PCM comparable to [5] 



and [10\ in Section 2.5 as the low-field limit of our general nonlinear theory. 



2.2. Cavity shape function s(r), and dependent energy A clLV 

Polarizable continuum models replace the liquid with a dielectric cavity surrounding 
the electronic system. In variants of the model suitable for treating solid surfaces 
(which typically require a plane- wave basis), the dielectric constant is smoothly 
switched from the vacuum value of 1 to the bulk liquid value ef, [HJ [9l [10] . The 

| Here and throughout this paper, we use atomic units Aire = e = h=m e =kg=l. 



Nonlinear fluid response in JDFT studies of battery systems 



5 



spatial modulation of the dielectric constant may be written as e(r) = 1 + (ej — l)s(r), 
so that s(r) £ [0, 1] describes the shape of the cavity. 

Further, these variants of PCM assume that the cavity shape s(r) = s(n(r)) is 
determined entirely by the local electron density. The exact functional form of s(n) 
is not important as long as it switches smoothly between at high electron densities 
and 1 at low electron densities, and rapidly approaches the extreme values away from 
the transition region. Following [6] , for the rest of this work, we use 



where the parameter n c sets the critical electron density around which the cavity 
smoothly 'switches on', and a controls the width of that transition. 

In the following subsections, we develop an ab initio theory for the nonlinear 
dielectric and ionic response of solvents, which we find to be the dominant effects 
at the charged or highly polar surfaces in electrochemical systems due to the strong 
electric fields. The cavity shape function, however, includes unknown parameters 
that are typically fit [7J [TO] to reproduce the solvation energies of small organic 
molecules. These solvation energies are sensitive to other free energy contributions 
such as cavitation and dispersion, which although negligible in the electrochemical 
systems of interest, cannot be ignored during the determination of fit parameters. 

These additional free energy contributions have complicated dependences on the 
shape of the cavity, for which several empirical approximations have been developed 
(see [7] for a review). Andreussi and coworkers |10) demonstrated that a simple 
empirical model expressing the sum of all these effects as an effective surface tension for 
the cavity works reasonably well for the solvation energies of small organic molecules. 
Since we need the additional effects only as auxiliary contributions during the fit to 
the molecular solvation data, we adopt their simplified model here by writing 



S 

where S is a surface area estimate for the cavity described by s(r), and r is an effective 
tension that is determined by the fit to solvation energies. 

2.3. Nonlinear dielectric internal energy, A £ 

The dielectric response of liquids includes contributions from molecular polarizability 
as well as rotations of molecules with permanent dipole moments. The response of 
highly polar solvents such as water is dominated by rotations. With increasing field 
strength, the molecular dipoles increasingly align with the electric field, eventually 
saturating the rotational response. The polarizability response, which includes 
electronic polarizability and flexing modes of the molecules, typically becomes stronger 
at higher fields. It is therefore important to consider all these contributions even for 
solvents whose linear response is dominated by rotations. 

The typical electric fields encountered in solvation can significantly saturate the 
rotational response of solvents, but are usually insufficient to access the nonlinear 
regime of the remaining contributions. We therefore split the internal energy of the 
dielectric A e into rotational A mt and polarization A po \ parts, and construct a nonlinear 
theory for the rotational part alone. 



s(n) = -erfc 
v ; 2 




ay/2 



(3) 




(4) 



Nonlinear fluid response in JDFT studies of battery systems 



G 



Within the polarizable continuum ansatz, the liquid consists of molecules 
distributed with the bulk density iV mo i modulated by the cavity shape function s(r). 
The internal energy corresponding to linear polarization response with an effective 
molecular polarizability Xmoi is 

A po i[P po i} = J driV mol s(r)^Xmoi-Ppoi( r ') I (5) 

where P po \(r) is the induced dipole moment per molecule. This dipole moment 
contributes a bound charge, p po i (r) = -V • (N mol s{r)P pol (r)). 

Physically, the nonlinearity of the rotational response arises from a competition 
between the rotational entropy of the molecules and their interaction with the self- 
consistent electric field. We therefore begin with the exact rotational entropy for an 
ideal gas of dipoles with the cavity-prescribed density N mo \s(r) at temperature T, 
then approximate rotational correlations, and write 



4ir Pe 



(6) 



A to t\pe,l] = J drTN mo is(r) J ^ Pe log Pe - l(r) f^J 

Here, p e (r) is the probability that a molecule at location r has its dipole oriented 
along unit vector e, and the second term of ^ constrains the normalization of p e (r) 
with Lagrange multiplier field l(r). 

The final term of ^ captures the correlations in dipole rotations within a 
'local polarization density approximation (LPDA)'. We choose the simplest possible 
form for this correction, quadratic in the local dimcnsionless polarization P rot (r) 
j j|p e (r)e, and constrain the prefactor a to reproduce the bulk linear dielectric 
constant e&. Finally, the rotational response contributes a bound charge p r ot( r ) = 
— V • ( Pmo \N mo is(r)P ro t(r)) within the local response approximation, where p mo i is 
the permanent molecular dipole moment. 

The Euler-Lagrange equation for minimizing the total free energy with respect 
to p e implies that, at equilibrium, the orientation distribution must be of the form 
Pe = cxp(e • e) for some vector field e(r). Using the remaining Euler-Lagrange 
equations to eliminate P po i(r) and Z(r) in favor of e(r), the sum of ^ and ([6| 
simplifies to 



A e [e(r)} - J drTN mol s(r) 



(7) 



with corresponding dielectric bound charge 

p £ (r) = -V • [ Pmo iN mol s(r)e (/(e) + X(l - a/(e)))] . (8) 

Here, /(e) = (e coth e — l)/e 2 is the effective dimensionless rotational susceptibility 
defined by P rot = /(e)e, and X = Xmoi^/Pmoi ^ s ^ s counterpart for the linear 
polarization response. 

The resulting theory for the dielectric has four solvent-dependent parameters 
(N mo i, Pmo i, X and a), of which the bulk molecular density, N mo \, is directly 
measurable. The effective molecule dipole moment in the liquid, p mo i, differs from the 
gas-phase value, and in principle, can be determined from measurements of the lowest 
order nonlinear response coefficient [19] . However, such measurements are difficult and 
not readily available for most solvents. Instead, we compute j? mo i as the self-consistent 
dipole moment of a single solvent molecule in a solvated ab initio calculation employing 
a nonlinear polarizable continuum description of the same solvent. The resulting dipole 
moment for the solvents used in this work are listed in table [3j Note that p mo i is larger 
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Figure 1. Comparison of the effective dielectric constant of water as a function of 
uniform externally applied field Erj for the nonlinear PCM model ^ compared to 
SPC/E molecular dynamics results 1211 and classical density functional predictions 
[15] . The effective dielectric constant is defined by e = Eq/E, where E = Eq—AttP 
is the net electric field including the screening due to the dielectric polarization 
density P. Within this theory, this response is determined entirely by bulk liquid 
properties e;,, eoo and iV mo i, along with the molecule dipole moment p mo i obtained 
from a self-consistent ab initio calculation solvated with the present model. 



than the gas phase dipole moment for all these solvents because charged centers in 
the molecule are surrounded by bound charges of the opposite sign which favor an 
increase in the polarization, as shown for water in figure [2j We also find that for 
water, p mo i = 0.94eao gratifyingly agrees with the SPC/E molecular dynamics model 
value of 0.92eao [20], in contrast to the gas phase value of 0.73earj. 

The remaining solvent-dependent parameters, the correlation factor for rotations 
a and the effective dimensionless polarizability X, are constrained to reproduce the 
bulk static and high frequency dielectric constants, th and eoo respectively. Using the 
bulk linear response of the above functional, we can analytically show that 

T(£oo - 1) 4?riV mo ip^ ol 
X — - — — = — and a = 3 — —. (9) 

since the rotational response freezes out and does not contribute to f m . In principle, 
eoo should be the dielectric constant at infrared frequencies between the rotational 
and vibrational resonances, but in the absence of experimental data in that frequency 
regime, we use the readily measurable optical dielectric constant, which is the square 
of the refractive index. 

We have therefore produced a density-functional theory for the nonlinear dielectric 
response of an arbitrary solvent constrained entirely by measurable macroscopic 
properties. The response at field strengths relevant for solvation is not accessible 
experimentally, but it has been estimated using molecular dynamics. Figure [l] 
demonstrates that the bulk nonlinear dielectric response of the present theory is in 
excellent agreement with molecular dynamics results for water [21] using the SPC/E 
pair potential model 20 . The present theory, which uses LPDA for rotational 
correlations, produces essentially the same nonlinear response to uniform electric 
fields as classical density functional theories with the scaled mean-field electrostatics 
approximation [22] . The minor differences between the present theory and the classical 
density functional results [15] shown in Figure [I] are due to electrostriction; the latter 
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theory accounts for changes in the equilibrium fluid density in the presence of a strong 
uniform electric field. 



2.4. Nonlinear ionic system internal energy, A K 

The previous section derived the dielectric response of liquids from the dipolar 
rotational and polarization response of liquid molecules to the local electric field. 
Ionic species in the liquid introduce Debye screening by contributing an additional 
monopolar response, which changes the local ionic density in response to the local 
electric potential. A simple description of this response at the linearized Poisson- 
Boltzmann level suffices for many electrochemical systems [5]. For the electrode- 
electrolyte interface, this level of theory corresponds roughly to the Gouy-Chapman- 
Stern model, but misses the nonlinear capacitance effects due to ion adsorption. 
Here, we explore whether a full Poisson-Boltzmann treatment within the polarizable 
continuum model ansatz captures these additional details. 

To represent the internal free energy of an ionic system comprising several species 
of charge Zi and bulk concentrations AT, each, we employ the exact expression for the 
ideal gas of point particle ions, approximating finite-size effects with a local density 
approximation, and write 



A K [{ m (r)}}=Tj2 J driV lS (r) 



jx(r) -x ) 2 
x (l - x ) 2 (l - x{r)Y ' 



fa(log V l) + l) + a(| (i)J(1 ((niJ 
Ideal gas 



Hard sphere 

The density of each ionic species is Nis(r)r)i(r), represented in terms of the 
enhancement rji(r) relative to the cavity prescription of ATjs(r). The charge-weighted 
sum of these densities contribute a net ionic bound charge p K = £\ ZiNis{r)rji(r). 



The first ideal gas term in ( 10 ) along with the mean-field electrostatic interaction 
in the third term of ([2| correspond to the Poisson-Boltzmann theory. That theory, 
however, does not limit the density of the ions in solution and presents an unphysical 
instability associated with infinite build-up of ions at regions of strong external 
potential. We resolve this instability by enforcing a packing limit on the ions via 



the second term of ( 10 ). This term captures local hard sphere correlations in terms of 
the packing fraction x(r) = ^ ViNiT]i(r), where V$ = 4nRf/3 is the volume per ion 
for each species (with ionic radius Ri). The functional form of this term is constrained 
to reproduce x(r) —> x^ = ^ ViNi in the bulk and to match the divergence in the 
equation of state of the hard sphere fluid 23; as x 1. 



2.5. Linear limit 

The free energy functional ([2| with dielectric free energy A € given by ^ and optional 
ionic free energy A K given by ( |10[ ) constitutes our nonlinear polarizable continuum 
model. Here, we show that the conventional linear polarizable continuum model is a 
limit of this more general theory. 

The rotational response from the dielectric is approximately linear when the 
energy scale of the molecular dipoles interacting with the field is much lower than 
the temperature (p mo i|V0| <C T), where 4>(r) is the total electrostatic potential. 
Similarly, the ionic response is approximately linear when Z\(f>\ -C T. Using the 
Euler-Lagrange equations to eliminate e(r) and /Lt(r) in favor of V0(r) and <j>{r) 
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respectively, expanding the free energy to quadratic order, and simplifying using ^ 
and the definition n 2 = 8?r N {on Z 2 /T, we find 



A e + A K = — I drs(r) 



\Vcp\ 2 2 cp 2 
[e& - 1)— h K 



(11) 



2 2 

with the corresponding total bound charge at linear order 

M»0 = ^ [(tb - 1) V • (s(r)V4>) - K 2 s(r)4>] . (12) 

The Euler-Lagrange equation for this simplified linear-response functional in 
terms of the single independent variable, 4>, can be rearranged into the familiar 
modified Poisson equation (or Hclmholtz equation for non-zero k) 

V 2 (j){r) + (e 6 - 1)V • (s(n(r))V0(r)) - k 2 s(n(r))4>(r) = -4np e i(r). (13) 

Finally, substituting the solution of ( 13 ) in the fluid free energy functional ^ with 
the dielectric and ionic energies given by (11 1, yields an equilibrium value for in 
the linear limit 

Air = 4~ + \ J dr Pcl (r) (V) - / dr'p^j . (14) 

Thus, the free energy functional approach to polarizable continuum models reduces, 
in the linear limit, to the standard approach [51 110) of replacing the vacuum Poisson 
equation with one modified by the fluid. 

2.6. Periodic systems and net charge 

An important class of applications of the nonlinear polarizable continuum model, and 
joint density-functional theory in general, is the study of electrochemical systems. 
These systems pose an interesting challenge as they often involve charged metal 
or doped semiconducting surfaces. The periodic boundary conditions necessary to 
accurately describe the delocalized electronic states of such systems complicate the 
addition of charge, since the energy per unit cell of a periodic system with net charge 
per unit cell is divergent. 

Including a counter electrode [Mj to keep the simulation cell neutral avoids 
this problem, but leads to wasted computational effort on irrelevant portions of the 
system and complicates the separation of physics at the two electrodes. Introducing 
Debye screening due to ions in the electrolyte neutralizes the unit cell with fluid 
bound charge and naturally captures the physics of the electrochemical double layer 
[S]. More importantly, unlike the Poisson equation obtained without ionic screening, 



the Hclmholtz equation for the electrostatic potential with screening ( 13 1 has a 
well-defined constant offset ('zero' of potential) in periodic boundary conditions. 
The resulting Kohn-Sham eigenvalues, and hence the electron chemical potential, 
correspond to a zero reference deep within the fluid, and this enables calibration of 
the electron chemical potential in DFT against electrochemical reference electrodes. 
(See [S] for details.) 

The electrostatic potential in the nonlinear polarizable continuum model is not 
obtained from a Helmholtz equation, and the bound charge in the ionic system does not 
neutralize a net charge in the electronic system at an arbitrary value of the independent 
variable /x(r). Here, we present the modifications required to correctly handle periodic 
systems within the nonlinear ionic screening model. 



Nonlinear fluid response in JDFT studies of battery systems 



10 



The mean-field Coulomb energy per unit cell of volume Q for the entire system 
with total charge density p tot = p e \ + p\ q can be written in the plane-wave basis as 
u = 2T2 Y,G^G\Ptot(G)\ 2 . Here, p to t{G) = J n drp tot (r) exp(~iG ■ r) for reciprocal 
lattice vectors G, and Kg = 47r/G 2 is the plane- wave basis Coulomb kernel. The 
divergent contribution at G = vanishes for neutral unit cells with Q to t = Ptot(0) = 0. 

The G~ 2 divergence results from the long-range 1/r tail of the Coulomb kernel. 
We can analyze the effect of the divergence by making the Coulomb kernel short- 
ranged on a length scale L much larger than the unit cell, and set L — > 00 at the end. 
The exact form of the regularization is not important; picking the Gaussian-screened 
potential crfc(r/L)/r results in the regularized Coulomb energy 

^=E^-t(G)| 2 + ^QL- (15) 

Note that the first term employs the standard Coulomb kernel in the plane-wave basis 
which drops the G — term by invoking a neutralizing background, and the second 
term contains the divergent part depending only on the total charge per unit cell. 

At finite large L, the second term of Ul penalizes Qtot 7^ and favors equilibrium 
configurations with small Qtot ■ The Euler-Lagrange equation for the net charge Qtot is 
A = dA/dQtot — —KL 2 Q tot /£l, where A is the total free energy excluding the divergent 
second term of (15 1. Note that dA/dQtot is finite for systems capable of adjusting 
their total charge, such as fluids with ionic screening, so that as L — > 00, Q to t 
in such a manner that A cx QtotL 2 remains finite. The absolute potential is also well 
defined in this situation with a G = contribution of 9C/i/9Qtot = 7rL 2 Q to t/^ = —A. 
Finally, note that 

U °° = E ^ol^°t(G)| 2 - AQtot (16) 

G^O 



results in the same Euler-Lagrange equation and equilibrium free energy as ( 15 ) in 
the L — > 00 limit, and therefore the divergent term in the Coulomb energy reduces to 
a charge-neutrality constraint imposed by Lagrange multiplier A. 

We incorporate this Lagrange multiplier constraint into the ionic free energy in 
plane-wave calculations, and retain the standard plane-wave Coulomb kernel with 
G = projected out for all electrostatic interactions. The constraint can be solved 
analytically for local nonlinear ions (section 2.4) in the commonly encountered case 
of a 'Z:Z' electrolyte consisting of two species of charge +Z and — Z (labeled with 
indices « = +,—) with bulk concentrations N{ on each. In this situation, we can show 
that substituting r]±(r) = exp(±(p + p±(r))), where p = —ZX/T is obtained by 
solving the neutrality constraint, reduces the constrained minimization over T]±(r) to 
an unconstrained minimization over ^±(1"). In particular, the neutrality constraint 
Q+e^° + Q-e-w + Q e] = yields 

Mo = log^^ — , (17) 

where Q± = ±N- mn Z J drs(r)e ± ^ 1 ^ and Q c \ = j drp c i(r) is the total charge of the 
electronic system. In this case, and for other joint density-functional theories which 
include ionic screening, the constraint contribution to <L4diei/<5p c i(?") in. the electron 
potential establishes the absolute reference for the Kohn-Sham eigenvalues and the 
electron chemical potential required for ah initio electrochemistry [S]. 
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2. 7. Implementation 

The nonlinear polarizable continuum model presented here and its linear counterpart 
have been implemented in the open source plane-wave electronic structure software 
JDFTx [35], designed for joint density- functional theory. The electronic density- 
functional theory segment of this software is based on conjugate gradients 
minimization [26J of an analytically continued total energy functional [27j , expressed in 
the DFT++ algebraic formulation [SB]. The fluid segment of JDFTx also employs the 
plane-wave basis and is discretized in the algebraic formulation for classical density- 
functional theories [15 ]. 

The valence electron density n(r) from standard pseudopotentials need to be 
augmented with a core electron density to prevent overlap of the fluid with the 
pseudopotential cores [5J. Hence, we compute the shape function using ([3| with 
^cavfr) — n(r) + ncoj-cfr), where n coro is the partial core density used for nonlinear 
core corrections [25] . 

The electrostatic interactions with the fluid involve the total charge density 
(both electronic and nuclear) of the material described in the electronic structure 
portion of the calculation, p e \{r) = n(r) + p nuc (r). Here, the nuclear charge density, 
p nuc (r) = — Zie ri ) l^ 2w ^/(27rui) 3 / 2 [§] is widened by a Gaussian resolvable on 
the charge density grid. The widened nuclear density is used only in the interaction 
with the fluid; the internal energies of the electronic system employ point nuclei in all 
terms. This width does not affect the interaction energy since the fluid and nuclear 
charge densities do not overlap, and the nuclear charge is spherically symmetric. 
However, it shifts the potential relative to the zero-width case, which we compensate 
exactly by adding the correction — 2irw 2 Zi/Cl to the electron potential, where O 
is the unit cell volume. 

Finally, regarding algorithms, the linear polarizable continuum models are 



minimized by solving the Helmholtz (or Poisson) equation ( 13 1 at every electronic 
iteration. Appropriate preconditioners for the involved linear conjugate gradients 
solver have been developed previously [5J . The free energy of the nonlinear polarizable 
continuum model is minimized using the Gummel iteration |30) . where the electronic 
system and the fluid are alternately minimized while holding the state of the other one 
frozen. This method is guaranteed to be globally convergent due to the variational 
principle, and typically converges adequately in 5-10 alternations for most systems 
studied. The fluid free energy A^iei is minimized with the scalar field (J,(r) and vector 
field e(r) as independent variables; the diagonal preconditioner in reciprocal space [jj] 



K^G) = 



Z(l-a/3)l 2 G 2 



Pmol 



(18) 



(G 2 + K 2 /e b ) 2 

for the p channel with the identity preconditioner on the e channel yields satisfactory 
convergence for the nonlinear conjugate gradients algorithm [26] . 



3. Results 



Strong electric fields at liquid interfaces typical of battery systems necessitate a theory 
for the nonlinear response of the liquid environment, such as the nonlinear polarizable 

§ Note that we employ an electron-is-positive charge convention, so that p nuc < and the charge 
of the electron is +1 in atomic units. 

|| This preconditioner is derived from an approximation to the Hessian of Adiel with respect to fJ.(r) 
and e(r) in the bulk linear limit. 
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continuum model of section [2] Section |3.2| calibrates the undetermined parameters of 
this theory against experimental solvation energies of molecules. For these molecules, 
and for metallic surfaces in section [373} we find results comparable to linear PCM's. 
However, for surfaces of ionic solids in section [3. 4| we find that inclusion of nonlinear 
effects are necessary in order to obtain qualitatively correct results. 

3.1. Computational Details 

We perform all calculations in this paper using the open source plane-wave density 
functional software JDFTx [25] at a plane wave cutoff of 30 Eh (1 Eh = 1 hartree 
w 27.21 eV). These calculations employ norm-conserving pseudopotentials generated 
by the Opium pseudopotential generator |31) with the PBE exchange and correlation 
functional [32 . The pseudopotentials for metal atoms include partial core corrections 
[2"5] , which are necessary to keep the fluid out of the pseudopotential cores as described 
in Section [2~7l 

The choice of exchange-correlation functional for molecular and surface systems 
is not straightforward |33) . and some argue that semi-local approximations can be 
inadequate for these systems [34j . Hybrid functionals which include exact exchange, 
or quantum Monte Carlo methods, are likely to be more accurate but are significantly 
more expensive than semi-local methods and hence unsuitable for rapid screening 
calculations. Here, we use the semi-local revTPSS exchange-correlation functional 
[35] which shows considerable promise for accurate calculations of surface phenomena 
including surface formation energies and molecular adsorption energies |36j . 

Molecular geometries for the calculations of section |3.2| are from the 
Computational Chemistry Comparison and Benchmark Database [37 . The surface 
geometries employed in sections |3.4| and |3.3| are constrained to the optimized bulk 
geometry for the central layer, while the remaining layers are fully relaxed for both 
the vacuum and fluid calculations. The fluid models assume ambient temperature 
T = 298 K for all calculations. 

3.2. Calibration to molecular solvation energies 

The nonlinear dielectric response of section |2.3| is completely constrained by ab initio 
and experimentally determined parameters, listed in table [3] for the solvents studied 
in this paper. However, the cavity shape function and the cavitation and dispersion 
terms, which are integral features of any polarizable continuum model, are unknown 
microscopic quantities that are typically constrained by a fit to solvation energies. 
Here, we fit the set of unknown cavity parameters for the nonlinear model and its 
linear limit to the same molecular solvation dataset using the same procedure, in 
order to facilitate a fair comparison between linear polarizable continuum models and 
our nonlinear theory. 

The molecular solvation dataset must contain experimental data that is both 
reliable and readily available. Organic molecules solvated in water satisfy this criterion 
and are commonly used in fitting parameters for polarizable continuum models [101 [7] . 
The molecules used in our fit are listed in figure [3J and the known solvent parameters 
for water are listed in table [3j Of the remaining parameters, we set the shape function 
width parameter a — 0.6 as in [6l [8] since the solvation energies are somewhat 
insensitive to it. We then determine the cavity transition electron density n c and the 
effective cavity tension r by a nonlinear least squares fit to the molecular solvation 
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n c (a d ) 


t (E h /4) 


RMS Error (kcal/mol) 


Nonlinear PCM 


1.0 x 10" 3 


9.5 x 10-° 


0.95 


Linear PCM 


3.7 x 10" 4 


5.4 x MT 6 


1.05 



Table 1. Fitted parameters of the nonlinear and linear polarizablc continuum 
models (PCM) and the corresponding RMS errors for solvation energies of the 
molecules listed in figure [3] 
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Figure 2. Bound charge in solvent water around a water molecule in the nonlinear 
and linear models. The smaller hydrogen atoms produce stronger fields on the 
solvent compared to the oxygen, resulting in much stronger saturation effects in 
the negative bound charge surrounding the hydrogens. In spite of the increased 
bound charge, the linear model yields approximately the same solvation energy 
as the nonlinear one due to compensation by the increased cavity size. 



energy dataset. 

The resulting fit parameters and optimized RMS error in solvation energy for the 
nonlinear and linear versions of the model are summarized in table [T] The smaller n c , 
and hence larger cavities, for the linear model as compared to the nonlinear one offset 
the overestimation of electrostatic interactions due to the lack of saturation effects. 
The lowered cavity tension r in the linear model then compensates for the increase in 
cavity area. Figure [2] demonstrates the consequences of these differences in the solvent 
bound charge surrounding a water molecule. The solvation energies predicted by the 
two models are in agreement as seen in figure [3j in spite of significantly larger bound 
charges in the linear case. Due to this cancellation, the linear model yields comparable 
accuracy to the nonlinear one for the solvation of organic molecules in water, but this 
is no longer the case when stronger electric fields come into play, as in some of the 
electrochemical systems we study next. 

3.3. Solvation of metallic surfaces 

Unlike the typical electrochemical interface, noble metal electrodes in electrolyte are 
less prone to complex chemical interactions at the surface, making them suitable 
candidates for an initial evaluation of our theory. Reactions are highly sensitive to 
the absolute electron chemical potential, which in experiments is typically reported 
relative to the standard hydrogen electrode (SHE). The absolute potential of the SHE 
relative to vacuum is difficult to establish experimentally; the estimates from different 




Figure 3. Solvation energies of molecules in water predicted by the nonlinear and 
linear polarizablc continuum models compared against the experimental values 
from 1551 1551 . 





VSHE (V) 


l^dip (V) 


RMS Error (V) 


Nonlinear PCM 


4.62 


0.46 


0.09 


Linear PCM 


4.68 


0.40 


0.09 



Table 2. Offset between theoretical and experimental PZC's, VshEi determined 
by a fit using the systems of figure |3J a), with corresponding RMS errors. Vshe 
represents the potential difference between an electron solvated deep in the fluid 
and the Standard Hydrogen Electrode. Vdip represents the potential due to 
the dipole moment at the fluid-metal interface, and is obtained as the difference 
between the theoretical PZC and the work function, averaged over the systems 
considered for each fluid model. 



experimental methods range from 4.4 V to 4.9 V 0D]. To make direct contact with 
experimental electrochemical observables, this experimentally uncertain quantity can 
be calibrated |S] in density-functional theory by comparing the theoretical chemical 
potentials for solvated neutral metal surfaces against the experimental potentials of 
zero charge (PZC's). The calibrations of the reference electrode potential within the 
linear and nonlinear models are remarkably similar, as shown in figure[4|a) and table[2] 
The absolute potential of zero charge includes contributions from the work 
function, which is essentially independent of the fluid theory, and from the dipole 
moment in the interfacial layers of the liquid. The minor differences in the calibrations 
of the two theories stem from this dipole moment contribution, as shown for aqueous 
electrolytes in table [2| The variation of surface charge with electrode potential is also 
similar for the two models, as shown for the solvated Pt(lll) surface of figure |2]jb) . 
In particular, the derivative of that variation, the so-called 'double-layer' capacitance, 
at the potential of zero charge is 14 and 15 /xF/cm 2 for linear and nonlinear PCM 
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Experimental PZC (V, relative to SHE) Electrode potential (V, relative to SHE) 

Figure 4. (a) Potentials of zero charge (PZC's) for the (111), (100) and (110) 
(left to right) each for silver (circles), gold (triangles), and copper (squares), 
predicted by the nonlinear and linear theories, compared to experiment |40| , The 
diagonal line for each theory compares theoretical and experimental values up to 
an overall fitted offset. (See table|2]) The silver and gold are solvated in aqueous 
1 M NaF electrolyte (ionic radii Na: 1.16 A F: 1.19 A), while the copper is in 
aqueous 1 M KCIO4 electrolyte (ionic radii K: 1.52 A and CIO4: 2.26 A), (b) 
Charge on a Pt(lll) surface in 1 M aqueous KCIO4 as a function of potential 
relative to the standard hydrogen electrode (SHE) for the two theories. 



respectively, which agrees well with an experimental estimate of 20 /iF/cm 2 [IT] for 
the above system. 

The agreement in the results of the linear and nonlinear theories demonstrated 
in figures |4]ja,b) and table [2] is due to the same cancellation of errors at play for 
solvation of molecules. The linear theory misses saturation in the rotational dielectric 
response, thereby overestimating it, yet compensates with an increase in cavity size. 
This cancellation of errors is possible since the typical magnitudes of electric fields 
under typical operating potentials are similar to those of the molecular shown 
in figure (7) 

Both models predict an approximately linear variation of surface charge with 
electrode potential (figure Qb)), which corresponds to a constant capacitance. This 
prediction contrasts with the experimental observation of a capacitance minimum at 
the potential of zero charge [UJ due to ion adsorption on the electrode surface. The 
formation of this so-called inner Helmholtz layer between the solid surface and the 
solvent is precluded by the cavity ansatz of polarizable continuum models. These 
details require either a higher level of theory capable of describing layering effects of 
ions such as a classical density-functional approach, or the inclusion of explicit ions 
into the quantum mechanical calculation. Nonetheless, both the linear and nonlinear 
PCM adequately describe the basic features of the ideal electrochemical interface, and 
are suitable for describing chemical reactions at metal electrode surfaces as long as all 
chemical bonds are treated quantum-mechanically. 

3.4- Solvation of ionic surfaces 

The surfaces of electrodes typically contain ionic compounds whose structure and 
composition vary with the chosen electrolyte and operating conditions. The details of 
this interface play a critical role in battery performance, and an accurate description 
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Solvent 


£6 




Pvac (eOo) 


Pmol (eoo) 


N mo i (% 6 


) 


t {E h /aD 


Water 


78.4 


1.78 


0.73 


0.94 


4.938 x 10- 


3 


9.5 x 10~ e 


DMC 


3.1 


1.87 


0.16 


0.16 


1.059 x 10- 


a 


2.05 x 10" 5 


THF 


7.6 


1.98 


0.69 


0.90 


1.100 x 10- 


3 


1.78 x 10" 5 


DMF 


38.0 


2.05 


1.50 


2.19 


1.153 x 10" 


3 


2.26 x 10~ 5 


PC 


64.0 


2.02 


1.97 


2.95 


1.039 x 10- 


3 


2.88 x 10~ 5 


EC 


90.5 


2.00 


1.93 


2.88 


1.339 x 10- 


3 


3.51 x 10" 5 



Table 3. Parameters describing water and commonly used lithium battery sol- 
vents, Dimethyl Carbonate (DMC), Tetrahydrofuran (THF), Dimethylformamide 
(DMF), Propylene Carbonate (PC) and Ethylene Carbonate (EC). The vac- 
uum dipoles (pvac) and self-consistent solvated dipoles p mo i are computed using 
density-functional theory as described in section |2.3| All remaining parameters 
are constrained by measured bulk properties 1441 . 




Figure 5. Solvation energies predicted by the nonlinear and linear models for 
surfaces of Li2 0, LiOH and LiF in the organic solvents from table [3] With 
increasing dielectric constant, the predictions of the linear model diverge from 
those of the nonlinear model due to missing saturation effects. This leads to 
qualitative differences here, unlike the case of the solvated molecules of figure [3] 
For some surfaces, the linear model suggests, perhaps incorrectly, that the flat 
surface is unstable by lowering the solvated surface energy relative to the bulk 
solid. 



of such surfaces in electrolyte environments is therefore crucial for modeling efforts 
towards improving battery systems. Reactions at the surface of a lithium metal anode, 
for example, can form Li 2 0, LiOH and LiF at the solid electrolyte interface [4121 143] . 
Here, we study these surfaces in contact with different organic solvents typical of 
battery systems as a testbed for fluid models applicable to battery systems. 

The solvents selected for this study are listed in table [3j Due to the dearth 
of experimental data for corresponding solvation energies, we here use the cavity 
shape parameters determined by the fit in section |3.2| We replace the effective 
tension r by the experimental surface tension, ignoring dispersion effects which are 
insignificant on the scale of the electrostatic energies in these highly polar systems. All 
remaining physical parameters that determine the dielectric response are constrained 



by experiment and ab initio calculations, as discussed in section 2.3 



The linear and nonlinear models predict similar solvation energies for the 
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Figure 6. Bound charge in solvent Ethylene Carbonate (EC) around a LiF (100) 
surface in the nonlinear and linear models, shown in a (011) slice. Saturation 
effects are stronger next to the smaller Li+ cations which produce significantly 
stronger fields on the solvent compared to the larger F — anions. In contrast to 
a water molecule solvated in water (Figure |2p, these effects are strong enough to 
qualitatively alter solvation energies, as shown in figure [5] 
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Figure 7. Effective rotational susceptibility at the average value of the 
dimensionless effective field e at the cavity surface (solvent-solute interface) for 
solvated molecules (circles), charged metal surfaces (triangles) and ionic surfaces 
(squares). The reduction in susceptibility due to saturation effects captured by 
the nonlinear model is missed by the linear one. Unlike the case of molecules and 
metal surfaces, the order of magnitude overestimation of the susceptibility by the 
linear model for ionic surfaces is not compensated by the increase in cavity size. 



aforementioned ionic compounds of lithium in solvents with low dielectric constants, as 
shown in figure [5] However, with increasing dielectric constant, the magnitude of the 
solvation energy increases more rapidly for the linear model, leading to disagreement 
by up to a factor of two for the most polar solvents. The linear model overestimates the 
electrostatic interaction due to a lack of saturation effects, but unlike in the molecular 
case, the increase in cavity size is insufficient to compensate for this error. In fact, 
for lithium fluoride in ethylene carbonate, as seen in figure [6j the linear model model 
overestimates the bound charge by an order of magnitude. Indeed, in this case, the 
result is a qualitative difference in the predicted stability of the solvated surface relative 
to the solid, with the linear model even predicting the solid to be thermodynamically 
unstable with respect to the formation of surfaces in this system. 
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The qualitative inadequacy of the linear model for ionic surfaces derives from 
the significantly stronger electric fields in these systems compared to solvated 
molecules and metallic surfaces. Figure [7] compares the average electric field and 
the corresponding rotational susceptibility at the solute-solvent interface for all 
the systems discussed above. The least polar neutral ionic surface still imposes 
a higher electric field than the most polar molecule or charged metallic surface 
at chemically relevant electrode potentials. The order of magnitude reduction in 
rotational susceptibility due to saturation effects in the ionic surfaces, compared to 
the modest reduction for the other systems, necessitates a nonlinear theory for the 
study of these types of battery systems. 

4. Conclusions 

Ab initio studies provide key insights into chemical processes in a wide range of 
systems, but have not yet approached battery chemistry with a realistic description 
of the electrolyte environment. Continuum solvation models provide an intuitive and 
computationally-efficient description of the environment and enable a focused study of 
the complex subsystems that require treatment at the electronic structure level. Our 
results indicate that standard polarizable continuum models fit to molecular solvation 
data perform poorly when applied to polar surfaces of the type often encountered at 
the SEI in battery systems. Consequently, one must exercise caution when attempting 
to apply standard solvation models available in both quantum chemical j7, 45, 46] and 
condensed matter [10l|47] ab initio software packages. As an alternative, the nonlinear 
theory presented here and implemented in |25j leverages the computational simplicity 
of the standard polarizable continuum models and extends their applicability to 
systems with the strong electric fields associated with ionic surfaces in electrochemical 
systems. 

The importance of nonlinear solvent response depends on the strength of electric 
fields at the interface, which in turn varies dramatically with system type, as 
highlighted in figure [7j For systems with moderate field strengths, such as the 
molecules and metal surfaces studied here, the linear models can compensate for 
the overestimated electrostatic response through an increase in cavity size. However, 
for systems with higher field strengths, such as ionic surfaces, this compensation is 
insufficient. The nonlinear polarizable continuum model developed here consistently 
describes all of these systems, and along with the technique developed in section |2.6| 
to determine the absolute electron chemical potential, enables electronic structure 
predictions for real electrochemical systems as a function of electrode potential. 

This theory provides a cost-effective yet accurate method for calculating 
properties of electrochemical systems of technological relevance, such as high- 
throughput screening potential battery materials. Cleaner surface experiments 
analogous to the solvation datasets available for molecules will help further refine the 
theory of solvation for these systems. In our study of battery electrode materials, we 
showcase our theory with electronic density functional calculations. The method can 
easily be used with other electronic theories such as coupled cluster or quantum Monte 
Carlo |14| . enabling the study of non-equilibrium properties. Using these techniques, 
we can include nonlinear solvation in transition state calculations important for 
understanding processes in energy systems, such as catalysis in fuel cells and ion 
diffusion on battery electrodes. 
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